The invention relates to signal processing methods and apparatus for performing convolution on data indicative of a pattern (e.g., image data indicative of a pixel array). In accordance with the invention, the convolution kernel is (or is approximated by) a spline function. Preferably the convolution kernel is (or is approximated by) a polynomial spline function that is piecewise a polynomial function.
Convolution is commonly performed on signals in many contexts, including the fields of sound, still image, video, lithography, radio (radar) signal processing. Typically, the signals to be convolved are pattern signals. Each of the expressions xe2x80x9cpatternxe2x80x9d and xe2x80x9cpattern signalxe2x80x9d is used herein in a broad sense to denote a one-dimensional sequence or two-dimensional (or higher dimensional) array of data words (which can be, but need not be pixels). Typically, the data words comprise binary bits, and the convolution is performed in discrete fashion on the binary bits using software, digital signal processing circuitry, custom hardware, or FPGA systems (field programmable gate array based computing systems).
The term xe2x80x9cdataxe2x80x9d herein denotes one or more signals indicative of data, and the expression xe2x80x9cdata wordxe2x80x9d herein denotes one or more signals indicative of a data word.
The motivations for implementing convolution rapidly, even when processing data indicative of very large patterns, are myriad. The present invention was motivated by the need for proximity correction in the field of lithography. In such problems, one attempts a two-dimensional convolution between data indicative of a large pattern xe2x80x9cpxe2x80x9d (where the pattern is a large one-dimensional or two-dimensional pixel array) and a diffusion kernel xe2x80x9cdxe2x80x9d. Often the kernel xe2x80x9cdxe2x80x9d is Gaussian or a superposition of Gaussians, or is otherwise a smooth kernel. More specifically, the present invention grew out of attempts to establish a suitable xe2x80x9cO(N)xe2x80x9d algorithm (an algorithm requiring not more than on the order of xe2x80x9cNxe2x80x9d multiplications and additions) for convolving a one-dimensional pattern comprising N pixels, where N is very large, with a Gaussian kernel (or other smooth kernel) such that the convolution is exact or very close to exact, or a suitable xe2x80x9cO(NNxe2x80x2)xe2x80x9d algorithm (an algorithm requiring not more than on the order of NNxe2x80x2 multiplications and additions) for convolving a two-dimensional pattern comprising NNxe2x80x2 pixels, where each of N and Nxe2x80x2 is very large) with a Gaussian kernel (or other smooth kernel) such that the convolution is exact or very close to exact.
The objective in performing proximity correction (in the field of lithography) is to generate a xe2x80x9crawxe2x80x9d optical signal (or xe2x80x9crawxe2x80x9d electron beam signal) which can be input to a set of reflective or refractive optics (or electron beam optics), in order to project the output of the optics as a desired pattern on a mask or wafer. To determine the characteristics of a raw optical signal (or raw electron beam signal) that are needed to cause projection of the desired pattern on the mask or wafer, a deconvolution operation is typically performed on a very large array of pixels (which determine a pattern xe2x80x9cpxe2x80x9d) in order to correct for the well known proximity problem. In the case of electron beam lithography, the proximity problem results from electron scattering in the substrate (mask or wafer) being written. Such scattering exposes broadened areas on the substrate to electrons (i.e., an area surrounding each pixel to be written in addition to the pixel itself), with the scattering effectively broadening the electron beam beyond the beam diameter with which the beam is incident on the substrate.
In nearly all proximity correction schemes, such a deconvolution operation includes at least one convolution step. Accordingly, in performing typical proximity correction, a very large array of pixels (determining a pattern xe2x80x9cpxe2x80x9d) must be convolved with a diffusion kernel. Although such a convolution is typically performed on a pattern comprising a very large array of binary pixels, this restriction is not essential in the following discussion and is not essential to implement the invention. Indeed, the invention can implement convolution on data indicative of any pattern xe2x80x9cpxe2x80x9d (which can be large or small, and can be either one-dimensional, two-dimensional, or more than two-dimensional, and can be determined by binary or non-binary data values), using any smooth convolution kernel xe2x80x9cdxe2x80x9d having characteristics to be described below.
For a data indicative of a pattern xe2x80x9cpxe2x80x9d and a convolution kernel xe2x80x9cdxe2x80x9d we consider the cyclic convolution:                     (                  d          xc3x97          p                )            n        =                  ∑                              i            +            j                    =                      n            ⁢                          xe2x80x83                        ⁢                          (                              mod                ⁢                                  xe2x80x83                                ⁢                N                            )                                          ⁢                        d          i                ⁢                  p          j                      ,
and an acyclic convolution which differs only in the indicial constraint and range:                     (                  d          ⁢                      xc3x97            A                    ⁢          p                )            n        =                  ∑                              i            +            j                    =          n                    ⁢                        d          i                ⁢                  p          j                      ,
where xe2x80x9cxAxe2x80x9d denotes that convolution operator x has an acyclic character.
In one-dimensional cases (in which the pattern p is an ordered set of N data values and the kernel is an ordered set of M values), the result of the cyclic convolution has length N (it comprises N data values), and the result of the acyclic convolution has length M+Nxe2x88x921. Both one- and two-dimensional scenarios are of interest. In the case of a two-dimensional pattern xe2x80x9cpxe2x80x9d (i.e., a pattern determined by an N by Nxe2x80x2 array of data values), the indices n, i, j and domain lengths N in the noted formulae are 2-vectors.
It is standard that the convolution dxc3x97p can be cast as an equivalent matrix-vector product:
dxc3x97pxe2x89xa1Dp,
where D is the circulant matrix of d (hereinafter the xe2x80x9ccirculantxe2x80x9d of d), whose 1-dimensional form is defined as:   D  =      (          xe2x80x83        ⁢                                        d            0                                                d            1                                                d            2                                                d            3                                    ⋯                                      d                          N              -              1                                                                        d                          N              -              1                                                            d            0                                                d            1                                                d            2                                    ⋯                                      d                          N              -              2                                                            ⋮                                      xe2x80x83                                                xe2x80x83                                                xe2x80x83                                                xe2x80x83                                    ⋮                                                  d            1                                                d            2                                                d            3                                                d            4                                    ⋯                                      d            0                                ⁢          xe2x80x83        )  
As will be shown later, conventional methods for convolution can be cast in the language of matrix algebra.
The central idea of the invention is to approximate a smooth kernel d by a polynomial spline kernel f (where f is a spline function f(x) which is piecewise a polynomial of degree xcex4 with M pieces fi(x)), and then to use appropriate operators that annihilate (or flatten) each polynomial of given degree (in a manner to be explained) to calculate the convolution of f and p quickly. In some embodiments, the smooth kernel d is approximated (in accordance with the invention) by a spline kernel f which is not a polynomial spline kernel, but which consists of M pieces defined over adjacent segments of its range (in typical two-dimensional cases, the latter spline kernel is a radially symmetric function whose range is some continuous or discrete set of values of the radial parameter). Though xe2x80x9csplinexe2x80x9d convolution in accordance with the invention has features reminiscent of conventional wavelet schemes and is an O(N) algorithm (as are wavelet schemes), an advantage of the inventive xe2x80x9csplinexe2x80x9d convolution is that it can be performed (on data indicative of a pattern p consisting of N data values) with cN arithmetic operations (multiplications and additions), whereas conventional wavelet convolution on the same data would require bN arithmetic operations, where the factor xe2x80x9cbxe2x80x9d is typically (i.e., with typical error analysis) significantly larger than the factor xe2x80x9cc.xe2x80x9d In other words, the implied big-O constant for the spline convolution is significantly smaller than the typical such constant for conventional wavelet convolution.
Other important advantages of the inventive spline convolution method will be discussed below. To better appreciate the advantages of the invention over conventional convolution methods, we next explain two types of conventional convolution methods: Fourier-based convolution and wavelet convolution.
As is well known, Fourier-based convolution is based on the elegant fact that if F is a Fourier matrix, say
Fjk=exe2x88x922xcfx80ijk/N,
then the transformation FDFxe2x88x921 of the circulant is diagonal, whence we compute:
Dp=Fxe2x88x921 (FDFxe2x88x921) Fp,
where the far-right operation Fp is the usual Fourier transform, the operation by the parenthetical part is (by virtue of diagonality) dyadic multiplication, and the final operation Fxe2x88x921 is the inverse Fourier transform. For arbitrary D one requires actually three Fourier transforms, because the creation of the diagonal matrix FDFxe2x88x921 requires one transform. However, if D is fixed, and transformed on a one-time basis, then subsequent convolutions Dp only require two transforms each, as is well known. The complexity then of Fourier-based cyclic convolution is thus O(N log N) operations (i.e., on the order of N log N multiplications and additions) for convolving a pattern p of length N (a pattern determined by N data values), for the 2 or 3 FFTs (Fast Fourier Transforms) required. It should be noted that the Fourier method is an exact method (up to roundoff errors depending on the FFT precision).
Another class of conventional convolution methods are wavelet convolution methods, which by their nature are generally inexact. The idea underlying such methods is elegant, and runs as follows in the matrix-algebraic paradigm. Assume that, given an N-by-N circulant D, it is possible to find a matrix W (this is typically a so-called compact wavelet transform) which has the properties:
(1) W is unitary (i.e. Wxe2x88x921 is the adjoint of W);
(2) W is sparse; and
(3) WDWxe2x88x921 is sparse,
where xe2x80x9csparsexe2x80x9d in the present context denotes simply that for sparse S, any matrix-vector product Sx, for arbitrary x, involves reduced complexity O(N), rather than say O(N2).
With the assumed properties, we can calculate:
Dp=Wxe2x88x921 (WDWxe2x88x921)Wp
by way of three sparse-matrix-vector multiplications, noting that unitarity implies the sparseness of Wxe2x88x921. Therefore the wavelet-based convolution complexity is O(N) for convolving a pattern p determined by N data values, except that it is generally impossible to find, for given circulant D, a matrix W that gives both sparsity properties rigorously. Typically, if the convolution kernel d is sufficiently smooth, then a wavelet operator W can be found such that within some acceptable approximation error the property (3) above holds. Above-noted properties (1) and (2) are common at least for the family of compact wavelets (it is property (3) that is usually approximate).
As noted, an advantage of xe2x80x9csplinexe2x80x9d convolution in accordance with the invention is that it can be performed (on data indicative of a pattern p comprising N data values) with cN arithmetic operations, whereas conventional wavelet convolution on the same data would require bN arithmetic operations, where (assuming typical error budgets) the factor xe2x80x9cbxe2x80x9d is significantly larger than the factor xe2x80x9cc.xe2x80x9d Among the other important advantages of the inventive xe2x80x9csplinexe2x80x9d convolution method (over conventional convolution methods) are the following: spline convolution is exact with respect to the spline kernel f, whereas wavelet convolution schemes are approximate by design. (and error analysis for wavelet convolution is difficult to implement); the signal lengths for signals to be convolved by spline convolution are unrestricted (i.e., they need not be powers of two as in some conventional methods, and indeed they need not have any special form); and spline convolution allows acyclic convolution without padding with zeroes (whereas conventional Fourier schemes require full zero padding for acyclic convolution).
In a class of embodiments, the invention is a method for performing cyclic or acyclic convolution of a pattern xe2x80x9cpxe2x80x9d (i.e., data indicative of a pattern xe2x80x9cpxe2x80x9d) with a smooth diffusion kernel d, to generate data indicative of the convolution result
r=Dp,
where D is the circulant matrix (sometimes referred to herein as the xe2x80x9ccirculantxe2x80x9d) of d. The pattern xe2x80x9cpxe2x80x9d can be one-dimensional in the sense that it is determined by a continuous one-dimensional range (or discrete set) of data values (e.g., pixels), or it can be two-dimensional in the sense that it is determined by a continuous two-dimensional range of data values (or array of discrete data values), or p can have dimension greater than two. In typical discrete implementations, the pattern p is one-dimensional in the sense that it is determined by a discrete, ordered set of data values (e.g., pixels) Pi, where i varies from 0 to Nxe2x88x921 (where N is the signal length), or it is two-dimensional in the sense that it is determined by an array of data values Pij, where i varies from 0 to Nxe2x88x921 and j varies from 0 to Nxe2x80x2xe2x88x921, or it has dimension greater than two (it is determined by a three- or higher-dimensional set of data values). Typically, the kernel d is determined by an array of data values dij, where i varies from 0 to Nxe2x88x921 and j varies from 0 to Nxe2x80x2xe2x88x921 (but the kernel d can alternatively be determined by a discrete set of data values d0 through dNxe2x88x921).
In preferred embodiments, the inventive method accomplishes the convolution Dp, by performing the steps of:
(a) approximating the kernel d by a polynomial spline kernel f (unless the kernel d is itself a polynomial spline kernel, in which case d=f and step (a) is omitted);
(b) calculating q=Bp=xcex94xcex4+1Fp, where F is the circulant of kernel f, and xcex94xcex4+1 is an annihilation operator (whose form generally depends on the degree xcex4 of the polynomial segments of F) which operates on circulant F in such a manner that xcex94xcex4+1F=B is sparse; and
(c) back-solving xcex94xcex4+1rapprox=q to determine rapprox=Fp.
In cases in which the kernel d is itself a polynomial spline kernel (so that d=f, and F=D), the method yields an exact result (rapprox=r=Dp)
Otherwise, the error inherent in the method is (fxe2x88x92d)xc3x97p, where xc3x97 denotes the cross product, and thus the error is bounded easily.
In one-dimensional cases (in which the pattern to be convolved is a one-dimensional pattern of length N), xcex94xcex4+1 has the form of the Nxc3x97N circulant matrix defined as follows:       Δ          δ      +      1        =      [          xe2x80x83        ⁢                                        +                          (                                                                                          δ                      +                      1                                                                                                            0                                                              )                                                            -                          (                                                                                          δ                      +                      1                                                                                                            1                                                              )                                                            +                          (                                                                                          δ                      +                      1                                                                                                            2                                                              )                                                            -                          (                                                                                          δ                      +                      1                                                                                                            3                                                              )                                                ⋯                          0                                      0                                      +                          (                                                                                          δ                      +                      1                                                                                                            0                                                              )                                                            -                          (                                                                                          δ                      +                      1                                                                                                            1                                                              )                                                            +                          (                                                                                          δ                      +                      1                                                                                                            2                                                              )                                                ⋯                          0                                      ⋮                                      xe2x80x83                                                xe2x80x83                                                xe2x80x83                                                xe2x80x83                                    ⋮                                                  -                          (                                                                                          δ                      +                      1                                                                                                            1                                                              )                                                            +                          (                                                                                          δ                      +                      1                                                                                                            2                                                              )                                                ⋯                                      xe2x80x83                                    0                                      +                          (                                                                                          δ                      +                      1                                                                                                            0                                                              )                                            ⁢          xe2x80x83        ]  
in which each entry is a binomial coefficient, and xcex4 is the maximum degree of the spline segments of spline kernel f. For example, xcex4=2 where the spline kernel f comprises quadratic segments. In two- or higher-dimensional cases, the annihilation operators can be defined as       Δ    =                            ∂                      x            1                                n            1                          ⁢                              ∂                          x              2                                      n              2                                ⁢                      xe2x80x83                    ⁢          …                    ⁢              xe2x80x83            ⁢              ∂                  x          d                          n          d                      ,
where ∂xnn is the n-th partial derivative with respect is to the h-th of d coordinates. For example, the Laplacian       ∇    2    ⁢      =                            ∂                      x            1                    2                ⁢                  xe2x80x83                ⁢        …            ⁢              xe2x80x83            ⁢              ∂                  x          d                          n          2                    
will annihilate piecewise-planar functions of d-dimensional arguments f=f (x1, x2, . . . xd). In the one-dimensional case, the end points of each segment (the xe2x80x9cpivot pointsxe2x80x9d) of spline kernel f may be consecutive elements di and di+1 of kernel d, and step (a) can be implemented by performing curve fitting to select each segment of the spline kernel as one which adequately matches a corresponding segment of the kernel d (provided that appropriate boundary conditions are satisfied at each pivot point, such as by derivative-matching or satisfying some other smoothness criterion at the pivot points).
In preferred implementations, step (c) includes a preliminary xe2x80x9cignitionxe2x80x9d step in which a small number of the lowest components of rapprox=Fp are computed by exact multiplication of p by a few rows of F, and then a step of determining the rest of the components of rapprox using a natural recurrence relation determined by the spline kernel and the operator xcex94xcex4+1. For example, in the one-dimensional case, the lowest components of rapprox are r0, r1, . . . , rxcex4, where xe2x80x9cxcex4xe2x80x9d is the maximum degree of the spline segments of spline kernel f (for example r0, r1, and r2 where the spline kernel comprises quadratic segments), and these xe2x80x9cxcex4xe2x80x9d components are determined by exact multiplication of p by xe2x80x9cxcex4xe2x80x9d rows of F. Then, the rest of the components xe2x80x9crxcex4xe2x80x9d are determined using a natural recurrence relation determined by the spline kernel and the operator xcex94xcex4+1. The xe2x80x9cignitionxe2x80x9d operation which generates the components r0, r1, . . . , rxcex4, can be accomplished with O(N) computations. The recurrence relation calculation can also be accomplished with O(N) computations.
In other embodiments, the inventive convolution method accomplishes the convolution r=Dp by performing the steps of:
(a) approximating the kernel d by a polynomial spline kernel f (unless the kernel d is itself a polynomial spline kernel, in which case d=f and step (a) is omitted);
(b) calculating q=Bp=xcex94xcex4Fp, where F is the circulant of kernel f and xcex94xcex4 is a flattening operator (whose form generally depends on the degree xcex4 of the polynomial segments of F, and which operates on circulant F such that B=xcex94xcex4F is almost everywhere a locally constant matrix); and
(c) back-solving. xcex94xcex4rapprox=q to determine rapprox=Fp. In one-dimensional cases (in which p has length N), xcex94xcex4 has the form of the Nxc3x97N circulant matrix:       Δ    δ    =            [              xe2x80x83            ⁢                                                  +                              (                                                                            δ                                                                                                  0                                                                      )                                                                        -                              (                                                                            δ                                                                                                  1                                                                      )                                                                        +                              (                                                                            δ                                                                                                  2                                                                      )                                                                        -                              (                                                                            δ                                                                                                  3                                                                      )                                                          ⋯                                0                                                0                                              +                              (                                                                            δ                                                                                                  0                                                                      )                                                                        -                              (                                                                            δ                                                                                                  1                                                                      )                                                                        +                              (                                                                            δ                                                                                                  2                                                                      )                                                          ⋯                                0                                                ⋮                                              xe2x80x83                                                          xe2x80x83                                                          xe2x80x83                                                          xe2x80x83                                            ⋮                                                              -                              (                                                                            δ                                                                                                  1                                                                      )                                                                        +                              (                                                                            δ                                                                                                  2                                                                      )                                                          ⋯                                              xe2x80x83                                            0                                              +                              (                                                                            δ                                                                                                  0                                                                      )                                                        ⁢              xe2x80x83            ]        -  
in which each entry is a binomial coefficient, and xcex4 is the maximum degree of the spline segments of spline kernel f. In higher-dimensional cases, the flattening operator xcex94xcex4 is defined similarly.
In the latter embodiments (those employing a flattening operator in step (b) rather than an annihilation operator), step (c) can be performed with fewer multiplications (since it can be performed with an xe2x80x9cignitionxe2x80x9d step which computes a smaller number of the lowest components of rapprox=Fp, i.e., only the components rapprox0, rapprox1, . . . , rapproxxcex4xe2x88x921, where xe2x80x9cxcex4xe2x80x9d is the maximum degree of the spline segments of spline kernel f) but it requires more addition operations.
In other embodiments, the inventive method accomplishes the convolution Dp (where D is the circulant of smooth kernel d) by performing the steps of:
(a) approximating the kernel d by a spline kernel f which is not a polynomial spline kernel (unless the kernel d is itself such a spline kernel, other than a polynomial spline kernel, in which case d=f and step (a) is omitted);
(b) calculating q=Bp=AFp, where F is the circulant of kernel f, and A is an annihilation or flattening operator, where A operates on circulant F in such a manner that AF=B is sparse when A is an annihilation operator, and A operates on circulant F in such a manner that AF=B is almost everywhere a locally constant matrix when A is a flattening operator); and
(c) back-solving Arapprox=q to determine rapprox=Fp.
In other embodiments, the invention is a computer programmed with software for performing convolution on data in accordance with any embodiment of the inventive method. Other embodiments of the invention include a digital signal processor including digital signal processing circuitry configured to perform convolution on data in accordance with any embodiment of the inventive method, an apparatus (such as custom or dedicated electronic circuitry, or a field programmable gate array based computing system (xe2x80x9cFPGA systemxe2x80x9d)) configured to perform convolution on data in accordance with any embodiment of the inventive method, and a lithography system including such digital signal processing circuitry, such custom or dedicated electronic circuitry, or such an FPGA system.
Also within the scope of the invention is a computer-readable storage medium which stores computer-executable instructions, wherein the instructions are such that a computer performs an embodiment of the inventive method in response to executing the instructions.